{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CD2QfkE-OV_q",
    "papermill": {
     "duration": 0.006583,
     "end_time": "2021-02-24T00:00:00.388859",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.382276",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Analyzing RCT with Precision by Adjusting for Baseline Covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 1367,
     "status": "ok",
     "timestamp": 1617484454603,
     "user": {
      "displayName": "Andy Haupt",
      "photoUrl": "",
      "userId": "04421689622138798294"
     },
     "user_tz": -120
    },
    "id": "pCj5gpZAOap8",
    "outputId": "ffd01f2c-a401-409d-ac17-258bcc10395b"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wCK49jmFOV_y",
    "papermill": {
     "duration": 0.005433,
     "end_time": "2021-02-24T00:00:00.399979",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.394546",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Jonathan Roth's DGP\n",
    "\n",
    "Here we set up a DGP with heterogenous effects. In this example, with is due to Jonathan Roth, we have\n",
    "$$\n",
    "E [Y(0) | Z] = - Z, \\quad E [Y(1) |Z] = Z, \\quad Z \\sim N(0,1).\n",
    "$$\n",
    "The CATE is\n",
    "$$\n",
    "E [Y(1) - Y(0) | Z ]= 2 Z.\n",
    "$$\n",
    "and the ATE is\n",
    "$$\n",
    "2 E Z = 0.\n",
    "$$\n",
    "\n",
    "We would like to estimate ATE as precisely as possible.\n",
    "\n",
    "An economic motivation for this example could be provided as follows: Let D be the treatment of going to college, and $Z$ academic skills.  Suppose that academic skills cause lower earnings Y(0) in jobs that don't require college degree, and cause higher earnings  Y(1) in jobs that require college degrees. This type of scenario is reflected in the DGP set-up above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2NSMUY6qOV_0",
    "papermill": {
     "duration": 0.182998,
     "end_time": "2021-02-24T00:00:00.589410",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.406412",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "n = 1000             # sample size\n",
    "Z = np.random.normal(size=n)         # generate Z\n",
    "Y0 = -Z + np.random.normal(size=n)   # conditional average baseline response is -Z\n",
    "Y1 = Z + np.random.normal(size=n)    # conditional average treatment effect is +Z\n",
    "D = np.random.binomial(1, .2, size=n)    # treatment indicator; only 20% get treated \n",
    "Y = Y1 * D + Y0 * (1 - D)  # observed Y\n",
    "Z = Z - Z.mean()       # demean Z\n",
    "data = pd.DataFrame({\"Y\": Y, \"D\": D, \"Z\": Z})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IEYPeB5tOV_1",
    "papermill": {
     "duration": 0.005616,
     "end_time": "2021-02-24T00:00:00.601042",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.595426",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Analyze the RCT data with Precision Adjustment\n",
    "\n",
    "Consider \n",
    "\n",
    "*  classical 2-sample approach, no adjustment (CL)\n",
    "*  classical linear regression adjustment (CRA)\n",
    "*  interactive regression adjusment (IRA)\n",
    "\n",
    "Carry out inference using robust inference, using the sandwich formulas (Eicker-Huber-White).  \n",
    "\n",
    "Observe that CRA delivers estimates that are less efficient than CL (pointed out by Freedman), whereas IRA delivers more efficient approach (pointed out by Lin). In order for CRA to be more efficient than CL, we need the CRA to be a correct model of the conditional expectation function of Y given D and X, which is not the case here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 515,
     "status": "ok",
     "timestamp": 1617484460691,
     "user": {
      "displayName": "Andy Haupt",
      "photoUrl": "",
      "userId": "04421689622138798294"
     },
     "user_tz": -120
    },
    "id": "xfV4PjpzOV_1",
    "outputId": "e4cd79bf-bae8-4812-ca70-83c4376a49da",
    "papermill": {
     "duration": 0.28228,
     "end_time": "2021-02-24T00:00:00.888966",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.606686",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "CL = smf.ols(\"Y ~ D\", data=data).fit()          \n",
    "CRA = smf.ols(\"Y ~ D + Z\", data=data).fit()      #classical\n",
    "IRA = smf.ols(\"Y ~ D + Z + Z*D\", data=data).fit() #interactive approach\n",
    "# we are interested in the coefficients on variable \"D\".\n",
    "print(CL.get_robustcov_results(cov_type=\"HC1\").summary())\n",
    "print(CRA.get_robustcov_results(cov_type=\"HC1\").summary())\n",
    "print(IRA.get_robustcov_results(cov_type=\"HC1\").summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2wHiWkMXOV_3",
    "papermill": {
     "duration": 0.007983,
     "end_time": "2021-02-24T00:00:00.905362",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.897379",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Using classical standard errors (non-robust) is misleading here.\n",
    "\n",
    "We don't teach non-robust standard errors in econometrics courses, but the default statistical inference for the `fit` procedure in python, `smf.ols()`, still uses 100 year old concepts, perhaps in part due to historical legacy.  \n",
    "\n",
    "Here the non-robust standard errors suggest that there is not much difference between the different approaches, contrary to the conclusions reached using the robust standard errors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 632,
     "status": "ok",
     "timestamp": 1617539918260,
     "user": {
      "displayName": "Andy Haupt",
      "photoUrl": "",
      "userId": "04421689622138798294"
     },
     "user_tz": -120
    },
    "id": "Ukoh8YM9OV_3",
    "outputId": "c20b6541-1de8-42ae-b9ee-8efeb5c4218c",
    "papermill": {
     "duration": 0.041385,
     "end_time": "2021-02-24T00:00:00.954907",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.913522",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(smf.ols(\"Y ~ D\", data).fit().summary())         \n",
    "print(smf.ols(\"Y ~ D + Z\", data).fit().summary())\n",
    "print(smf.ols(\"Y ~ D + Z + Z*D\", data).fit().summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hoT263l7OV_4",
    "papermill": {
     "duration": 0.010058,
     "end_time": "2021-02-24T00:00:00.974929",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.964871",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Verify Asymptotic Approximations Hold in Finite-Sample Simulation Experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_execution_state": "idle",
    "_uuid": "051d70d956493feee0c6d64651c6a088724dca2a",
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 19916,
     "status": "ok",
     "timestamp": 1615754191853,
     "user": {
      "displayName": "Andy Haupt",
      "photoUrl": "",
      "userId": "04421689622138798294"
     },
     "user_tz": -60
    },
    "id": "QQZG7z_kOV_4",
    "outputId": "e3d0a6fc-2133-4f1d-ca75-33160d59246f",
    "papermill": {
     "duration": 3.55396,
     "end_time": "2021-02-24T00:00:04.538599",
     "exception": false,
     "start_time": "2021-02-24T00:00:00.984639",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "\n",
    "np.random.seed(123)\n",
    "\n",
    "def exp(it, n):\n",
    "    np.random.seed(it)\n",
    "    Z = np.random.normal(size=n)\n",
    "    Y0 = -Z + np.random.normal(size=n)\n",
    "    Y1 =  Z + np.random.normal(size=n)\n",
    "    D = np.random.binomial(1, .2, size=n)\n",
    "    Y = Y1 * D + Y0 * (1-D)\n",
    "\n",
    "    Z = Z - Z.mean()\n",
    "    data = pd.DataFrame({\"Z\": Z, \"D\": D, \"Y\": Y})\n",
    "    CL = smf.ols(\"Y ~ D\", data).fit().params[\"D\"]\n",
    "    CRA = smf.ols(\"Y ~ D + Z\", data).fit().params[\"D\"]\n",
    "    IRA = smf.ols(\"Y ~ D + Z+ Z*D\", data).fit().params[\"D\"]\n",
    "    return CL, CRA, IRA\n",
    "\n",
    "n = 1000\n",
    "B = 1000\n",
    "res = Parallel(n_jobs=-1, verbose=3)(delayed(exp)(it, n) for it in range(B))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "mLfkZGwCdBWj"
   },
   "outputs": [],
   "source": [
    "res = np.array(res)\n",
    "CLs, CRAs, IRAs = res[:, 0], res[:, 1], res[:, 2]\n",
    "print(\"Standard deviations for estimators\")\n",
    "print(np.sqrt((CLs**2).mean()))\n",
    "print(np.sqrt((CRAs**2).mean()))\n",
    "print(np.sqrt((IRAs**2).mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "DONE_analyzing-rct-with-precision.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 7.390855,
   "end_time": "2021-02-24T00:00:04.661485",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-02-23T23:59:57.270630",
   "version": "2.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
